Skewness Test (D’Agostino)#

The skewness test checks whether a sample shows statistically significant asymmetry.

It is commonly used as a diagnostic for:

  • symmetry assumptions (many parametric procedures implicitly assume roughly symmetric errors)

  • whether a transformation (log/sqrt/Box–Cox) might be useful

  • as one component of normality diagnostics (but it is not a full normality test by itself)

Learning goals#

By the end you should be able to:

  • explain what skewness measures (right/left tail intuition)

  • state the hypotheses behind the skewness test

  • implement D’Agostino’s skewness test with NumPy only

  • interpret the test statistic and p-value (including one-sided alternatives)

  • understand common pitfalls (outliers, sample size, “significant but tiny” effects)

import numpy as np
import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots

pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")

np.set_printoptions(precision=6, suppress=True)

rng = np.random.default_rng(42)

Prerequisites#

  • Sample mean and central moments

  • Z-scores and p-values

  • The standard normal distribution

If any of these are fuzzy, you can still follow the notebook: the code is written to make each step explicit.

1) Intuition: what skewness measures#

Skewness is a signed measure of asymmetry:

  • Positive skewness → a longer right tail (“right-skewed”)

  • Negative skewness → a longer left tail (“left-skewed”)

A common (moment-based) sample skewness is

[ ext{skew}(x) = rac{m_3}{m_2^{3/2}}, \qquad m_k = rac{1}{n}\sum_{i=1}^n (x_i - ar{x})^k. ]

This uses the third central moment (m_3), so it is very sensitive to outliers.

def sample_skewness(x: np.ndarray) -> float:
    "Moment-based sample skewness: m3 / m2^(3/2) using 1/n moments."
    x = np.asarray(x, dtype=float).ravel()
    x = x[np.isfinite(x)]
    n = x.size
    if n < 3:
        return float("nan")

    mean = x.mean()
    centered = x - mean
    m2 = np.mean(centered**2)
    if m2 == 0:
        return 0.0
    m3 = np.mean(centered**3)
    return float(m3 / (m2 ** 1.5))


n_demo = 4000
x_sym = rng.normal(0, 1, size=n_demo)
x_right = rng.exponential(scale=1.0, size=n_demo)
x_left = -rng.exponential(scale=1.0, size=n_demo)

skews = {
    "Symmetric (normal)": sample_skewness(x_sym),
    "Right-skewed (exponential)": sample_skewness(x_right),
    "Left-skewed (-exponential)": sample_skewness(x_left),
}

skews
{'Symmetric (normal)': 0.003983194961544123,
 'Right-skewed (exponential)': 1.979806008441165,
 'Left-skewed (-exponential)': -1.901827954731168}
fig = make_subplots(
    rows=1,
    cols=3,
    subplot_titles=[f"{k}<br>skew={v:.3f}" for k, v in skews.items()],
)

datasets = [
    ("Symmetric (normal)", x_sym),
    ("Right-skewed (exponential)", x_right),
    ("Left-skewed (-exponential)", x_left),
]

for col, (name, x) in enumerate(datasets, start=1):
    fig.add_trace(
        go.Histogram(
            x=x,
            nbinsx=80,
            histnorm="probability density",
            name=name,
            marker=dict(line=dict(width=0)),
        ),
        row=1,
        col=col,
    )

fig.update_layout(
    height=320,
    width=1000,
    showlegend=False,
    title_text="Skewness: symmetric vs right/left skew",
)
fig.update_xaxes(range=[-5, 5])
fig.show()

2) What the skewness test is (and what it is not)#

Goal#

We want to know whether the skewness we observe could plausibly be due to sampling noise when the underlying distribution is normal (which has population skewness 0).

Hypotheses#

For the standard two-sided skewness test:

  • (H_0): the population skewness is the same as the normal distribution (i.e. 0)

  • (H_1): the population skewness is not 0

You can also run one-sided versions:

  • alternative='greater': (H_1): skewness (> 0) (right-skew)

  • alternative='less': (H_1): skewness (< 0) (left-skew)

What it is not#

  • It is not a complete normality test. A distribution can be symmetric (skew ≈ 0) but still very non-normal (e.g. heavy-tailed).

  • It does not tell you why data are skewed (mixtures, censoring, outliers, bounded support, etc.). You should always visualize the data.

3) D’Agostino’s skewness test: the core idea#

D’Agostino’s skewness test takes the sample skewness (b_2) and transforms it into a value (Z) that is approximately standard normal under (H_0):

[ Z pprox \mathcal{N}(0,1) \quad ext{(when the data are normal and } n\ge 8 ext{)}. ]

That lets us compute a p-value just like a z-test.

The transformation (as used by SciPy’s skewtest)#

Let (n) be the sample size and (b_2) the (moment-based) sample skewness. Define

[ y = b_2\sqrt{ rac{(n+1)(n+3)}{6(n-2)}}. ]

Then compute

[ eta_2 = rac{3,(n^2 + 27n - 70)(n+1)(n+3)}{(n-2)(n+5)(n+7)(n+9)}, \quad W^2 = -1 + \sqrt{2(eta_2 - 1)}. ]

and

[ \delta = rac{1}{\sqrt{ frac{1}{2}\ln(W^2)}}, \quad lpha = \sqrt{ rac{2}{W^2 - 1}}. ]

Finally,

[ Z = \delta,\operatorname{asinh}!\left( rac{y}{lpha} ight) = \delta,\ln\left( rac{y}{lpha} + \sqrt{\left( rac{y}{lpha} ight)^2 + 1} ight). ]

Why the sample size constraint?#

The approximation is derived under normality and is intended for (n\ge 8). For smaller samples, you should prefer resampling-based checks or simply treat skewness as an effect size.

4) NumPy-only implementation#

The functions below implement:

  • sample skewness (b_2)

  • D’Agostino’s transformed statistic (Z)

  • p-values using a NumPy approximation to the normal CDF (via an erf approximation)

The goal here is transparency, not squeezing out every last bit of numerical accuracy.

def erf_approx(x: np.ndarray) -> np.ndarray:
    "Vectorized erf approximation (Abramowitz & Stegun 7.1.26)."
    x = np.asarray(x, dtype=float)

    sign = np.sign(x)
    ax = np.abs(x)

    p = 0.3275911
    t = 1.0 / (1.0 + p * ax)

    a1 = 0.254829592
    a2 = -0.284496736
    a3 = 1.421413741
    a4 = -1.453152027
    a5 = 1.061405429

    poly = (((((a5 * t + a4) * t + a3) * t + a2) * t + a1) * t)
    y = 1.0 - poly * np.exp(-ax * ax)
    return sign * y


def norm_cdf(z: np.ndarray) -> np.ndarray:
    "Standard normal CDF Φ(z) using erf approximation."
    z = np.asarray(z, dtype=float)
    return 0.5 * (1.0 + erf_approx(z / np.sqrt(2.0)))


def pvalue_from_z(z: np.ndarray, alternative: str = "two-sided") -> np.ndarray:
    "Convert a z-statistic to a p-value under N(0,1)."
    z = np.asarray(z, dtype=float)

    if alternative == "two-sided":
        p = 2.0 * norm_cdf(-np.abs(z))
    elif alternative == "greater":
        p = norm_cdf(-z)
    elif alternative == "less":
        p = norm_cdf(z)
    else:
        raise ValueError("alternative must be one of: 'two-sided', 'greater', 'less'")

    return np.clip(p, 0.0, 1.0)


def dagostino_z_from_skewness(b2: np.ndarray, n: int) -> np.ndarray:
    "D'Agostino transform from sample skewness b2 to Z (approx N(0,1) under H0)."
    if n < 8:
        raise ValueError("D'Agostino skewness test requires n >= 8")

    b2 = np.asarray(b2, dtype=float)
    nf = float(n)

    y = b2 * np.sqrt(((nf + 1.0) * (nf + 3.0)) / (6.0 * (nf - 2.0)))

    beta2 = (
        3.0
        * (nf**2 + 27.0 * nf - 70.0)
        * (nf + 1.0)
        * (nf + 3.0)
        / ((nf - 2.0) * (nf + 5.0) * (nf + 7.0) * (nf + 9.0))
    )

    w2 = -1.0 + np.sqrt(2.0 * (beta2 - 1.0))
    delta = 1.0 / np.sqrt(0.5 * np.log(w2))
    alpha = np.sqrt(2.0 / (w2 - 1.0))

    # asinh(u) == log(u + sqrt(u^2 + 1)) but is more numerically stable
    return delta * np.arcsinh(y / alpha)


def dagostino_skewtest(x: np.ndarray, alternative: str = "two-sided") -> dict:
    "D'Agostino skewness test (NumPy-only implementation)."
    x = np.asarray(x, dtype=float).ravel()
    x = x[np.isfinite(x)]
    n = x.size
    if n < 8:
        raise ValueError("D'Agostino skewness test requires at least 8 observations")

    b2 = sample_skewness(x)
    z = dagostino_z_from_skewness(b2, n)
    p = pvalue_from_z(z, alternative=alternative)

    return {
        "n": int(n),
        "skewness": float(b2),
        "z": float(np.asarray(z)),
        "pvalue": float(np.asarray(p)),
        "alternative": alternative,
    }

5) Worked examples + interpretation#

Interpretation guide (two-sided):

  • Small p-value (e.g. < 0.05) → evidence the data are asymmetric (skewness ≠ 0) relative to normal sampling noise.

  • The sign of Z tells you the direction:

    • Z > 0 → right-skew (long right tail)

    • Z < 0 → left-skew (long left tail)

Remember: with large (n), even tiny skewness can become “statistically significant”. Always look at the effect size (the skewness value) and the plot.

def summarize_test(x: np.ndarray, name: str, alternative: str = "two-sided") -> dict:
    res = dagostino_skewtest(x, alternative=alternative)
    return {
        "dataset": name,
        "n": res["n"],
        "skewness": res["skewness"],
        "z": res["z"],
        "pvalue": res["pvalue"],
        "alternative": res["alternative"],
    }


n = 60
x1 = rng.normal(0, 1, size=n)
x2 = rng.exponential(1.0, size=n)
x3 = -rng.exponential(1.0, size=n)

results = [
    summarize_test(x1, "Normal sample"),
    summarize_test(x2, "Right-skew (exponential)"),
    summarize_test(x3, "Left-skew (-exponential)"),
]

results
[{'dataset': 'Normal sample',
  'n': 60,
  'skewness': 0.33403897478489747,
  'z': 1.1391470464804974,
  'pvalue': 0.25464193619392184,
  'alternative': 'two-sided'},
 {'dataset': 'Right-skew (exponential)',
  'n': 60,
  'skewness': 1.0681651963774472,
  'z': 3.204306467065742,
  'pvalue': 0.001354010697767749,
  'alternative': 'two-sided'},
 {'dataset': 'Left-skew (-exponential)',
  'n': 60,
  'skewness': -2.314018629424382,
  'z': -5.399094653879995,
  'pvalue': 6.713003297686981e-08,
  'alternative': 'two-sided'}]
fig = make_subplots(
    rows=1,
    cols=3,
    subplot_titles=[
        f"{r['dataset']}<br>skew={r['skewness']:.3f}, Z={r['z']:.2f}, p={r['pvalue']:.3g}" for r in results
    ],
)

for col, (x, r) in enumerate([(x1, results[0]), (x2, results[1]), (x3, results[2])], start=1):
    fig.add_trace(
        go.Histogram(
            x=x,
            nbinsx=30,
            histnorm="probability density",
            name=r["dataset"],
            marker=dict(line=dict(width=0)),
        ),
        row=1,
        col=col,
    )

fig.update_layout(
    height=320,
    width=1100,
    showlegend=False,
    title_text="Skewness test on small samples (n=60)",
)
fig.show()

One-sided alternatives#

If you have a directional question, you can use a one-sided test:

  • greater: “is it right-skewed?” (skewness > 0)

  • less: “is it left-skewed?” (skewness < 0)

One-sided p-values are smaller only when the observed skew is in the hypothesized direction.

x = rng.exponential(1.0, size=60)

res_two = dagostino_skewtest(x, alternative="two-sided")
res_greater = dagostino_skewtest(x, alternative="greater")
res_less = dagostino_skewtest(x, alternative="less")

{
    "two-sided": res_two,
    "greater": res_greater,
    "less": res_less,
}
{'two-sided': {'n': 60,
  'skewness': 2.616686399904669,
  'z': 5.78336467838644,
  'pvalue': 7.344590824409636e-09,
  'alternative': 'two-sided'},
 'greater': {'n': 60,
  'skewness': 2.616686399904669,
  'z': 5.78336467838644,
  'pvalue': 3.672295412204818e-09,
  'alternative': 'greater'},
 'less': {'n': 60,
  'skewness': 2.616686399904669,
  'z': 5.78336467838644,
  'pvalue': 0.9999999963277046,
  'alternative': 'less'}}

6) Why sample size matters (p-value vs skewness)#

Because the test statistic is a standardized measure of skewness, larger samples make it easier to detect small asymmetries.

The plot below shows the two-sided p-value as a function of sample skewness (b_2) for different (n).

b2_grid = np.linspace(-2.5, 2.5, 501)

fig = go.Figure()

for n in [10, 30, 100, 500]:
    z_grid = dagostino_z_from_skewness(b2_grid, n=n)
    p_grid = pvalue_from_z(z_grid, alternative="two-sided")
    fig.add_trace(go.Scatter(x=b2_grid, y=p_grid, mode="lines", name=f"n={n}"))

fig.add_hline(y=0.05, line_dash="dash", line_color="black")

fig.update_layout(
    title="Two-sided p-value vs sample skewness (D'Agostino transform)",
    xaxis_title="sample skewness b2",
    yaxis_title="p-value",
    yaxis_type="log",
    height=420,
    width=900,
)
fig.show()

7) Does Z really look standard normal under H0?#

Under the null (normal data), the transformed statistic (Z) is designed to be approximately (\mathcal{N}(0,1)).

We can sanity-check that with a small Monte Carlo simulation.

def skewness_vectorized(samples: np.ndarray) -> np.ndarray:
    "Vectorized skewness for samples shaped (n_sims, n)."
    samples = np.asarray(samples, dtype=float)
    mean = samples.mean(axis=1, keepdims=True)
    centered = samples - mean
    m2 = np.mean(centered**2, axis=1)
    m3 = np.mean(centered**3, axis=1)

    out = np.zeros_like(m2)
    mask = m2 > 0
    out[mask] = m3[mask] / (m2[mask] ** 1.5)
    out[~mask] = 0.0
    return out


n = 20
n_sims = 10_000
samples = rng.normal(0, 1, size=(n_sims, n))

b2 = skewness_vectorized(samples)
z = dagostino_z_from_skewness(b2, n=n)

z_mean = float(np.mean(z))
z_std = float(np.std(z, ddof=0))
rejection_rate = float(np.mean(pvalue_from_z(z, alternative="two-sided") < 0.05))

{"z_mean": z_mean, "z_std": z_std, "rejection_rate@0.05": rejection_rate}
{'z_mean': 0.010093399552716806,
 'z_std': 0.9999562884320329,
 'rejection_rate@0.05': 0.0476}
# Histogram of Z with a standard normal PDF overlay

x_grid = np.linspace(-4, 4, 400)
pdf = (1.0 / np.sqrt(2 * np.pi)) * np.exp(-0.5 * x_grid**2)

fig = go.Figure()
fig.add_trace(
    go.Histogram(
        x=z,
        nbinsx=60,
        histnorm="probability density",
        name="Simulated Z (H0: normal)",
        opacity=0.75,
        marker=dict(line=dict(width=0)),
    )
)
fig.add_trace(go.Scatter(x=x_grid, y=pdf, mode="lines", name="N(0,1) PDF", line=dict(color="black")))

fig.update_layout(
    title=f"Under H0, Z is approximately N(0,1) (n={n}, sims={n_sims})",
    xaxis_title="Z",
    yaxis_title="density",
    height=420,
    width=900,
)
fig.show()

8) Practical notes and pitfalls#

  • Outliers can dominate: skewness uses ((x-ar{x})^3). One extreme point can flip the conclusion.

  • Statistical vs practical significance: with large (n), tiny skewness can yield tiny p-values.

  • Not a full normality test: skewness ≈ 0 does not imply normality.

  • Independence matters: if observations are dependent (time series, clustered data), p-values can be misleading.

  • Use visuals: combine the test with histograms/QQ-plots and domain knowledge.

9) Optional: compare with SciPy#

If you have SciPy available, you can verify that the NumPy implementation matches scipy.stats.skewtest.

x = rng.lognormal(mean=0.0, sigma=0.8, size=200)

ours = dagostino_skewtest(x)

try:
    from scipy.stats import skewtest

    scipy_res = skewtest(x)
    comparison = {
        "ours": ours,
        "scipy": {"z": float(scipy_res.statistic), "pvalue": float(scipy_res.pvalue)},
        "abs_diff_z": float(abs(ours["z"] - scipy_res.statistic)),
        "abs_diff_p": float(abs(ours["pvalue"] - scipy_res.pvalue)),
    }
except Exception as e:
    comparison = {"ours": ours, "scipy": None, "error": repr(e)}

comparison
{'ours': {'n': 200,
  'skewness': 2.5556762139133857,
  'z': 9.304601238442498,
  'pvalue': 0.0,
  'alternative': 'two-sided'},
 'scipy': {'z': 9.304601238442498, 'pvalue': 1.3449606122843268e-20},
 'abs_diff_z': 0.0,
 'abs_diff_p': 1.3449606122843268e-20}

Exercises#

  1. Simulate data from a symmetric but heavy-tailed distribution (e.g. Student t) and see how often the skewness test rejects.

  2. Create a mixture of two normals with different means and observe how skewness changes.

  3. For a fixed skewness, increase (n) and see how quickly the p-value shrinks.

References#

  • R. B. D’Agostino et al. (1990), A suggestion for using powerful and informative tests of normality

  • SciPy documentation for scipy.stats.skewtest